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ABSTRACT 

For a decade, N-body simulations have revealed a nearly universal dark matter density profile, 
which appears to be robust to changes in the overall density of the universe and the underlying power 
spectrum. Despite its un iversality, the physic al origin of this profile has not yet been well understood. 
Semi-analytic models bv lBarnes et al.l (|2005l ) have suggested that the density structure of dark matter 
halos is determined by the onset of the radial orbit instability (ROI). We have tested this hypothesis 
using N-body simulations of collapsing dark matter halos with a variety of initial conditions. For 
dynamically cold initial conditions, the resulting halo structures are triaxial in shape, due to the mild 
aspect of the instability. We examine how variations in initial velocity dispersion affect the onset of 
the instability, and find that an isotropic velocity dispersion can suppress the ROI entirely, while a 
purely radial dispersion does not. The quantity a 2 jv\ is a criterion for instability, where regions with 
<r 2 /v 2 < 1 become triaxial due to the ROI or other perturbations. We also find that the radial orbit 
instability sets a scale length at which the velocity dispersion changes rapidly from isotropic to radially 
anisotropic. This scale length is proportional to the radius at which the density profile changes shape, 
as is the case in the semi-analytic models; however, the coefficient of proportionality is different by a 
factor of ^2.5. We conclude that the radial orbit instability is likely to be a key physical mechanism 
responsible for the nearly universal profiles of simulated dark matter halos. 

Subject headings: dark matter — galaxies: halos — galaxies:formation — galaxies: evolution — insta- 
bilities 
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1. INTRODUCTION 



The universal density profile of dark matter ha- 
los has been w idely observe d in cosmological N- 

body simulations (ICole fe Lacevfl 996; M oore et aljl999t 
Bullock et~aT1 1200H Ivan den Bosch! 120021: iNavarro et all 
2004f k however, the physical origins are not yet under- 
stood. Whether formed in an isolated collapse or by hier- 
archical merging, dark matter halo density profiles have 
a characteristic double-power law shape. This shape is 
usually discussed in terms of the slope of the density pro- 
file 7 = d(logp)/c?(logr). The canonical halo profile has 
7 ~ — 1 in the inner regions and 7 ~ —3 in the outer 
regions jP ubinski & Carlbcrg 1991; Nav arro et al.lll996l 
119971 iHuss et all 119991: iMacMillan et~all I2006T ). While 
there are disagreements on the exact values of 7, the 
halo density profiles appear to be similar over decades 
in radius. There is also evidence that such profiles may 
accurat ely fit the luminous parts of early-type galax ies 
as well (|Dalcanton fc HogaiJl200ll : iMerritt et aHl2005D . 

Within dark matter halos, quantities other 
than density show universal profile s as well. 
For example. iTavlor fc Navarrol (|2001h and later 
iDehnen fc McLaughlin! ( 20051 ) showed that the phase- 
space density, represented by p/er 3 , has a constant 
power-law slope of a ~ 1.9 for over two and a half 
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decades in radius in s pite o f a varying density profile 
(see also I Austin et al.l (|2005l ) for an analytic exploration 
of this behavior). There is also evidence for a univer- 
sality in the relation between the slope of the density 
profile and the velocity anisotropy. The anisotropy 
parameter, is defined as 



= 1- a 2 /2a 2 



(1) 



where er^ is the velocity dispersion in the tangential di- 
rection and a r is the velocity dispersion in the radial 
direction. describes the degree of anisotropy of the 
velocity dispersion in a spherical halo, which in simula- 
tions ranges from isotropic in the core (0 = 0) to radi- 
ally anisotropic at the outskirts (0 = l) (Cole fc Lacevl 
19961: ICarlberg et alJ Il997t iFukushige fc Makinol l2001h 



Hansen fc Moord (|2006h find a direct correlation between 



and the density profile slope for halos with a wide vari- 
ety of initial conditi ons (i.e. isolated collapse, mergers). 
lHansen et al.l (|2006f ) find another universality, this time 
in the velocity distribution function, for a similar variety 
of initial conditions. 

The radial orbit instability (ROI) may be the cause 
of the apparent universality of these dark matter halo 
properties. The ROI occurs in anisotropic spherical sys- 
tems composed of particles with predominantly radial 
orbits. The instability arises when particles in process- 
ing elongated loop orbits experience a torque due to 
a slight asymmetry. This torque causes them to lose 
some angular momentum and move towards the cen- 
ter of the system. Particles with small enough angular 
momenta become trapped in box orbits, aligning them- 
selves with the growing bar and reinforcing the initial 
perturbation. As a result, the halo becomes triaxial 
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in sha pe. This phenomeno n has been examined in de- 
tail bv iDeionghe fc Merrittl (|1988| ) using N-body simu- 
lations, bvlWeinberd (119911) using linear analysi s, and by 
IHuss et al.l (|1999l) and iMacMillan et all (|2006l ) for the 
specific case of isolated dark matter halo collapse. These 
groups find that the onset of the ROI corresponds to a 
flattening of the central density cusp of the halo, and 
may be responsible for the double power-law shape of 
the halo. The ROI is reviewed in Merritt (1999; §6.2). 

Semi-analytic models have also been used to ex- 
amine the sphe rical collapse of dar k matter halos, 
begi nning with iGunn fc Gottl (119721). These m od- 
els (|Goti) [l975t iGmTnl fl977t iBertschingerl fl985l ) in- 
clude only radial motions and produce single power- 
law density profiles. When non-radial motions are 
included, halo density profiles range fr om power-laws 
(iRvden fc Gunnl [198% to NF W-like (iHiotelid [200l 
He Delliou fc Henriksenl l2003t lAscasibar et alj l2004f) . 
Recently, however. IBarnes et al.l (120051) ex tended the 
semi-analytic model of Williams ct al. (2004) to include 
a physical representation of the ROI. The resulting halos 
have a double power-la w density slop e simil ar to those 
of N-body simulations. IBarnes et alJ (|2005f ) concluded 
that the density scale length and the anisotropy radius 
are correlated, and that the ROI is directl y responsible 
for th e shape of the density profile. While IBarnes et all 
(2005:) 's work is suggestive, the complexities of a dynami- 
cal instability are not easily captured by analytic or semi- 
analytic techniques. We therefore turn to studying the 
link between the ROI and halo structure by using high 
resolution N-body simulations. 

In this paper, we examine in detail the effect of random 
motions on the onset of the ROI and the final structure 
of collapsed dark matter halos, and attempt to verify 
the relation betwee n the scale length an d the anisotropy 
radius reported bv IBarnes et alJ (|2005h . We analyze a 
set of N-body simulations of isolated, collapsing dark 
matter halos with a variety of initial velocity dispersions 
to study the evolution of halo properties. Considering a 
range of velocity dispersions allows us to supress the ROI 
in some cases, helping us to isolate its physical effects. 
In §2 we describe our simulations and show that they are 
robust to resolution and softening effects. We describe 
the properties of the halos in §3, including the shape and 
anisotropy evolution, the effects of velocity dispersion, 
and the scale-length - anisotropy radius relation. We 
summarize our results in §4. 

2. SIMULATIONS 



We performed simulatio ns using PKDGRAV (|Stadell 
120011: IWadslev et all 120041 ) . a parallel KD Tree gravity 
solver. The initial system is an isolated, spherical halo 
with a Gaussian radial density profile, mimicking an over- 
density in the early universe. All particles are given Hub- 
ble flow velocities, such that the halo is initially expand- 
ing. We evolve the system in physical coordinates for 
12.8 Gyr, which was set to allow 83% of the mass of the 
halo to collapse by the present time; this corresponds to 
a starting redshift of ~ 12 accordin g to the most recent 
results from WMAP (|Spergel et al 



r g to t r 
. 2007). 



We ran several simulations with comparable density 
profiles but with a range of initial velocity dispersions, 
allowing us to estimate the threshold at which the ROI 
will occur. In principle, a dynamically "warmer" system 



should resist the ROI (|Merritt fc Aguilarl Il985l ). since 
any perturbations will be washed out by the tangential 
velocity dispersions. We assigned initial random veloci- 
ties (ct) to the particles assuming a Gaussian distribution 
with a mean of zero; these numbers are then added to the 
existing Hubble flow velocities. The amplitudes of these 
initial velocity dispersions were la, 2a, and 3a, where a 
is 18% of the circular velocity at the virial radius of the 
final system. Simulations with a higher velocity disper- 
sion are too warm to collapse at all. We also carried out 
simulations with no velocity dispersion, at both standard 
and high resolution. Finally, we attempted to isolate the 
importance of tangential velocity dispersions by carry- 
ing out a simulation with a pure radial initial velocity 
dispersion of 3a. Table 1 lists the initial value of the 
ratio — 2T/W (see M3 . 3[) for each simulation for a more 
intuitive grasp of the value of each a. 

In addition to various velocity dispersions, we ran sim- 
ulations with a range of softenings and particle resolu- 
tions. Our highest resolution simulation has over 2.8 
million particles, while the simulations testing the effects 
of velocity dispersion have about one fifth of that num- 
ber. These standard-resolutio n simulations are compa - 
rable in number of particles to IMacMillan et al] (T2 006). 
but ar e much larger than the simulations in IHuss et all 
(1999), which have ^10,000 particles within the virial 
radius of each halo. The softening, e, was set to 0.054 
times r2oo of the final collapsed halo, and is the same 
for all simulations. Time steps were set to an accuracy 
criterion proportional to rjy/e/a where rj = 0.2 is the 
timestep criterion an d a is the acceleration of a particle 
(|Wadslev et al.ll200l . We use a force accuracy criterion 
of 0.55. A summary of all simulations is in Table 1. 

2.1. Resolution and Softening Tests 

Figure Q] shows how the structure of our simulated ha- 
los are affected by particle number (also referred to in 
the text as resolution) . The top panels show the density 
profile for the fully evolved, Oct halo at standard reso- 
lution (left) and at high resolution (right). Each pro- 
file is fit with a NFW profile (purple dotted line) and a 
iNavarro et alJ (|2004D profile (also known as an Einasto 
profile) (red dashed line). The profiles for the two runs 
are similar over a large range in radius. However, as 
can be seen by the residuals from the fits, the stan- 
dard resolution simulation does not adequately model 
the inner core of the halo (r < 0.05r2oo, where we define 
r2oo as the virial radius in which the enclosed density 
is 200 times the critical density of the universe. There 
is also a difference in the innermost cores when looking 
at the z — axis ratio profiles (Fig. [U middle panels). 
The high-resolution run (right) shows a somewhat more 
triaxial core, particularly at the center. However, the 
anisotropy and phase-space density profiles show no sig- 
nificant differences between the two simulations (Fig. [IJ 
bottom panels, left and right, respectively). We ana- 
lyze these profiles in more detail in §3. The bulk of 
the difference between the two resolutions is apparent 
only within ~ 0.05 x r2oo- Outside of this radius our 
medium-resolution simulations are robust to resolution 
effects, while inside this radius the profiles may be less 
certain, and biased towards shallow density profiles and 
rounder shapes. 

We have also tested the impact of gravitational soften- 
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TABLE 1 
Simulation Properties 



Simulation 


# particles 


Velocity Dispersion 3 


Softcning b 


Initial -2T/W 


V0.H 


2865813 





1 


1.44 


V0.M 


573302 





1 


1.45 


V1.M 


573302 


1 


1 


1.49 


V2.M 


573302 


2 


1 


1.59 


V3.M 


573302 


3 


1 


1.64 


V3r_M 


573302 


3 (radial only) 


1 


1.87 


V0.M.S1 


573302 





2 


1.45 


V0.M.S2 


573302 





0.5 


1.45 



a in units of cr, where a is 18% of the circular velocity of the "cold' 
"cold" standard-resolution halo at z = 0. 

ing on our simulations. We performed two simulations, 
identical to the standard-resolution, Oct halo but with 
softening increased/decreased by a factor of two. De- 
creasing the softening length does not noticeably change 
any important quantities (density, velocity dispersion, 
axis ratios, phase-space density). Increasing the soft- 
ening length results in a shallower core density profile, 
due to the decreased gravitational force at small inter- 
particle distances. The phase-space density profile is 
similarly affected by the softening changes; however, the 
anisotropy profile shows no difference when the soften- 
ing lengths are changed. These effects are only noticeable 
within a radius of 0.05r2oo, the same radius within which 
we may be affected by resolution limitations (0.054r2oo)- 
We conclude that our results are robust beyond a radius 
of 0.05r2oo- This distance corresponds to 10 kpc for a 
Milky Way-size halo. 

3. RESULTS 

As the halo collapses, we measure the density profile, 
velocity dispersion tensor, and the axis ratios as a func- 
tion of radius at each time step. We expect to see higher 
central densities and larger radial velocity dispersions 
with increasing time. However, if the ROI is present we 
also expect to see deviations from spherical symmetry 
develop in regions with large radial velocity dispersions. 
These deviations would be accompanied by increasing 
tangential velocities due to the additional torques. 

3.1. Axis Ratios 

The most obvious evidence for the ROI is the growth 
of triaxiality in the shape of the halo. We see this be- 
havior clearly in our simulations. We measure axis ra- 
tios as a function of radius by dividing the halo into 196 
equal particle-number radial bins based on particle den- 
sity, and then c alculating the mom e nt-of- inertia tensor 
for each bin (see iBabul fc Starkmanl (|1992t ) for an appli- 
cation of this method). The square root of the ratios of 
the principal moments of inertia in the directions of the 
three axes of symmetry give the axis ratios of the parti- 
cles within each radial bin. This method more accurately 
tracks the radial variations of the axis ratios than meth- 
ods which analyze all particles interior to a given radius. 
This latter method is always dominated by the shape of 
the inner halo, where the majority of the virialized halo's 
mass is located, and thus is insensitive to variations in 
the axis ratios at large radii. 

After 0.6 Gyr, the first shell 5 reaches the center lead- 

5 While N— body halos do not have specific "shells" like those in 



standard-resolution halo at z = 0. b in units of 0.054 times T200 °f the 

ing to a large radial velocity dispersion. In response, a 
triaxial structure develops in the core of the halo. The 
resulting axis ratio profiles can be seen in Figure[21 which 
portrays the evolution of the axis ratios b/a and c/a at a 
range of timesteps for the high-resolution, zero velocity 
dispersion halo. The center of the halo quickly becomes 
triaxial and remains so throughout the collapse. The tri- 
axiality is fairly prolate (c ~ 0.726) and appears as a cen- 
tral bar. The outer edges of the halo have yet to collapse, 
and so remain spherical. As the collapse progresses, the 
triaxial region of the halo grows outwards, and eventu- 
ally the majority of the halo settles to an equilibrium b/a 
ratio of ~ 0.6. 

If the formation of the triaxial bar in Figure [2] is due 
to the ROI, then it should be suppressed when the or- 
bits are less radial. Our simulations with larger isotropic 
initial velocity dispersions would then be expected to 
show weaker triaxiality. The particles in these simula- 
tions have larger tangential velocities, which increases 
the typical angular momentum of particles' orbits, sup- 
pressing the formation of box orbits, and thus decreas- 
ing the strength of the ROI. Figure [3] shows the result- 
ing fully evolved axis ratio profiles for a series of sim- 
ulations with increasing initial random velocities. The 
standard and high resolution runs with zero velocity dis- 
persion show the central triaxial signature of the ROI. 
The la halo also shows some of this structure, but with 
a rounder core; this halo appears to undergo a weak form 
of the ROI. The 2a and 3a halos, on the other hand, re- 
main nearly spherical throughout until the end of their 
evolution. Their final shapes are characterized by rather 
spherical cores, with the inner axis ratios becoming pro- 
gressively more spherical with increasing initial tangen- 
tial velocities. A pre-existing isotropic velocity distribu- 
tion therefore prevents the onset of self-gravitating in- 
stabilities such as the R OI. Our results are supported by 
iTrenti fc Bertinl (|2006f ). who find that an isotropic core 
acts to stabilize a radially collapsing system against the 
ROI, even when the halo is highly anisotropic overall. 
Although we do not show the full time evolution, the ha- 
los that develop a strong or moderate ROI do so quickly, 
establishing a prolate triaxial core within 1 Gyr of the 
onset of collapse. Thus the triaxial shape indicates the 
presence of an instability; a phenomenon such as random 
scattering of the particles off of the potential would not 
create any preferred shape. 

We note that the increase in radial velocity dispersion 

the semi-analytic models, we use this terminology here to refer to 
particles at a similar initial radius and gravitational potential. 
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Fig. 1. — Comparison of the standard— and high— resolution fully evolved halos. Top row: Density profiles of the standard resolution run 
(left) and the high resolution run (right). Fits are NFW with a concentration of 14.6 (purple dotted line) and the Einasto profile with a 
= 0.2 (red dashed line). The residuals depict the superior core— fitting in the high-resolution simulation and also show that the Einasto 
profile is a good fit to r ~ 10r200- Middle row: Axis ratio profiles for the standard run (left) and the high-res run (right). Black line is b/a, 
and red dashed line is c/a. Bottom left: Anisotropy profiles for the standard (red dashed) and high-res (black solid) simulations. Bottom 
right: phase-space density profiles for the standard (red dashed) and high-res (black solid) runs. The standard resolution run shows slight 
deviations in density profile at r < 0.05r2oo an d is slightly less triaxial. However, the anisotropy and phase-space density profiles are 
similar at high and standard resolutions. 



alone has no significant effect on the final shape of the 
halo compared to the Oct collapse. In the bottom right 
panel of Figure [H we plot results for a simulation with 
large initial radial velocity dispersions, but no tangen- 
tial velocities. Physically, the increased initial velocity 
dispersion does little to change the final radial velocity 
dispersion of the core, since the latter is dominated by 
the larger radial motions produced by the gravitational 
collapse of the halo as a whole. Thus, the resulting axis 
ratio profile is nearly identical to the case with no initial 
velocity dispersion. The only significant difference is a 
slight delay in the collapse of the outer regions. 

To determine the physical likelihood of having velocity 
dispersions as high as that in the 3ct halo, we compared 
our N-body simulati ons to the semi-analytic model of 
IWilliams et al.l (|2004h . which treats an evolving pertur- 
bation as a series of gravitationally interacting shells. 



The model includes the effects of nearby substructure, 
which imparts a velocity dispersion onto the collapsing 
shells. This velocity dispersion increases with time, as 
the local substructure becomes more developed and a 
shell's time of exposure to the substructure increases. In 
contrast, in our N-body simulations of isolated collaps- 
ing halos the velocity dispersions are assigned only at the 
outset. 

The semi-analytic model shows that a velocity disper- 
sion of 2ct is given to a shell detaching from the Hubble 
flow at a redshift of 6.7, while 3ct is assigned to a shell at 
z ~ 4.3. These redshifts correspond to times when the 
ROI is occurring or has just occurred in the Oct N-body 
simulations; at these times the innermost particles' mo- 
tions are dominated by the gravitational potential of the 
halo and will not be strongly affected by external torques 
due to substructure. Subsequent infalling particles will 
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0.01 0.10 1.00 10.00 100.000.01 0.10 1.00 10.00 100.00 

r/r 200 r/r 200 

Fig. 2. — Axis ratio profiles of the highest resolution halo at various times. Black solid lines are the ratio b/a, red dashed lines are the 
ratio c/a. The halo profiles quickly become prolate in the inner regions, with the prolate region expanding to larger radii with increasing 
time. The profiles are normalized to r2oo a t each timcstcp. We interpret the various "bumps" in the profiles as being crossing shells that 
have not yet reached equilibrium. 



have increased velocity dispersion due to these effects of 
substructure; however, the ROI will have already pro- 
duced a triaxial core at this point. 

To summarize, the particles which collapse first and 
undergo the ROI may be slightly affected by torques due 
to substructure, but these forces are not strong enough 
to cause the large velocity dispersions which quench the 
ROI. In addition, our 2ct and 3ct halos which do not un- 
dergo the ROI are not physically plausible, since the lo- 
cal substructure cannot impart such large torques on the 
innermost particles at the early times at which they col- 
lapse. This suggests that in a cosmological context, the 
ROI will have time to operate in the central regions, be- 
fore cosmological torques have a significant effect. 

3.2. Velocity Dispersion Evolution 

If the ROI is responsible for the triaxiality in the inner 
regions of low velocity dispersion halos (Figure 13), then 



we should expect the onset and radial extent of the triax- 
iality to be associated with regions of highly anisotropic 
radial velocity dispersions and growing tangential veloci- 
ties. For our simulations, we calculate oy and at each 
timcstcp using 196 spherical bins. Velocity dispersions 
are determined by subtracting the average bulk velocity 
from the velocities in each bin. 

In Figures 2] and [5] we show the radial variation of the 
radial and tangential velocity dispersions as a function 
of time for our different halos. The central radial veloc- 
ity dispersion of the Oct halo undergoes a sharp increase 
as the innermost regions collapse into the center of the 
halo (Figure QJ upper right panel, black solid line) . The 
radial velocity dispersion increases outward with time 
as the halo collapses and material from larger radii falls 
inward and virializes. The Ict halo (red dotted line) un- 
dergoes a slower collapse and thus the jump in ct> is de- 
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0.01 0.10 1.00 10.00 100.000.01 0.10 1.00 10.00 100.00 



200 J-' 1 200 

Fig. 3. — The final axis ratio profiles of each simulation. Black solid lines are the ratio b/a, red dashed lines are the ratio c/a. As the 
velocity dispersion is increased, the halo becomes increasingly spherical in the inner regions (r < 0.2r2oo)i suppressing the ROI. 



layed. The slower collapse is due to the larger dynamical 
pressure support provided by the additional initial veloc- 
ity dispersion. However, despite these initial variations 
in collapse times, after a few billion years every halo has 
collapsed and virialized, and all of the ay profiles become 
nearly indistinguishable. This likeness is not surprising 
considering the similarities in the final density and other 
profiles (see §3.4, 3.5). 

The tangential velocity dispersion also shows an in- 
crease in the center of each halo early in the evolution 
(Figured}. In the case of the Oct halo, this increase is due 
to the torques applied by the radial orbit instability. In 
the halos which do not undergo the ROI, ct^ still increases 
due to the increased mass in the center of the halo, as 
the cores grow from material infalling from large radii. 
This increase in the central mass will cause any orbiting 
particles to migrate inwards due to the larger gravita- 
tional forces. However, as the particles move inwards, 
their angular velocity must increase to conserve angular 



momentum. Thus, halos with high initial ct^ can also in- 
crease their tangential velocity dispersion, even when the 
ROI does not operate. This effect will also contribute to 
increasing ct^ in the Oct halo once the ROI has imparted 
net angular momentum to some of the particles . As 
with the radial velocity dispersion, the final tangential 
velocity dispersions are very similar for each halo, even 
though the mechanisms to create such profiles arc differ- 
ent, and the resulting halos have different shapes. The 
similarity suggests that after the ROI initially increases 
ct (p , all subsequent evolution is set by angular momen- 
tum conservation during collapse. Therefore, while the 
ROI does not uniquely create the observed universal halo 
properties such as the NFW density profile, it does play 
a crucial role in seeding the initial tangential velocities 
in the dynamically "cold" collapse case, which then set 
the subsequent evolution of the halo. 

Armed with the velocity dispersion evolution shown in 
Figures 2] and we now return to the origin of the tri- 
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Fig. 4. — The evolution of the radial velocity dispersion profile for each of the halos. Simulations represented are no dispersion (black 
solid line), la dispersion (red dotted line), 2<r dispersion (green long— dashed line), and the 3<r dispersion (purple dashed line). The timestep 
is labelled on each panel. Initially, oy is a small fraction of the final virial velocity. As the system virializes, however, a r /v c ,200 tends to 
approach unity. The virialization clearly occurs from the inside out. 



axiality seen for all halos in Figure [3] (i.e. in the inner 
regions of the Oct and Ict halos, and the outer regions 
of the 2(7 and 3ct halos). In Figure [6] we plot the ratio 
of the square of the velocity dispersion over the circular 
velocity at a given radius. When the velocity dispersions 
are less than the circular velocity (i.e. a 2 jv 2 < 1), ha- 
los become unstable to triaxial perturbations. If a slight 
perturbation occurs, the tangential velocity dispersion is 
not sufficient to wash out the perturbation in less than 
an orbital time. This is the case at early times for the Oct 
and Ict halos in the inner regions, as well as in the outer 
regions for the 2a and 3ct halos. Where this mismatch 
of velocity scales occurs, any triaxial perturbation will 
grow until the ratio of a 2 jv 2 approaches unity. Com- 
paring Figure [6] to Figure [3] shows that indeed the outer 
prolateness seen in the high-cr halos begins at the radius 
where a 2 /v 2 drops below 1. Interior to this point, the 
velocity dispersions are high enough that any triaxiality 



is smoothed out by random motions before it can be re- 
inforced. In the Iow-ct halos, a 2 /v 2 starts well below 1, 
but the ROI acts quickly to raise the velocity dispersions 
of the centers of the halos, resulting in triaxial cores. 

3.3. Global Anisotropy 

The detailed behavior of the velocity dispersions in 
§3.2 can be simplified to a single measure of the velocity 
structure. A commonly used parameter for assessing the 
velocity structure of a halo is the global anisotropy: 

2T r /T =< a 2 > I < a 2 > (2) 

where T r and are kinetic energies of particles in the 
radial and tangential directions, respectively, and the av- 
erages are performed over the entire system. This quan- 
tity will be heavily weighted by the structure of the in- 
ner halo, which contains most of the halo mass. The 
global anisotropy parameter is intended to be used to 
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timestep is labelled on each panel. 



describe spherical equilibrium systems, and we note that 
the most complete way to compute the anisotropy of 
an evolving, nonsphcrical halo is to make a complete 
three dimensional map of the velocity dispersion ten- 
sor. This quantity is difficult to represent, and we find 
that when we calculate the global anisotropy parameter 
it closely follows the central velocity dispersion tensor. 
Thus we use the quantity 2T r /T^ to characterize the typ- 
ical anisotropy content at each timestep, even though our 
halos are evolving and deviate from spherical symmetry. 
Increasing values of the global anisotropy parameter in- 
dicate an increasing prevalence of radial orbits, and for 
a spherical equilibrium system, an increased likelihood 
of undergoing the ROI. Though a specific threshold for 
onset of the ROI has not been agreed upon, and depends 
moderately on the initial distribution function. Studies 
have found the thr eshold to be as lo w as TF r lT& ~ 1.4 
(iBarnes et al.lll986f) or as high as 2.5 (|Merritt fc Aguilarl 
119851: iMeza h Zamoranolll997l ) for equilbrium systems. 



At each timestep, the values < of > and < er^ > for 
the entire halo are found by averaging the squares of the 
values of ay and 0$ (see ^3.2[) from every radial bin nor- 
malized by the number of particles in each bin. The re- 
sulting values of 2T r /T r j ) vs. time are presented in the up- 
per panel of Figure for the various initial velocity dis- 
persions. At 0.5 Gyr, most halos have 2T r /T^ < 0.5 (the 
Oct and 3cr radial-only halos have negligibly small tangen- 
tial motions at this time, causing the value of 2T r /T$ to 
be extremely large). However, the subsequent evolution 
of the profiles is quite different. The low velocity halos 
that undergo the ROI show strong increases in 2T r /T^, 
which peaks at 0.6 Gyr, when the large infall velocities 
produce large values of < of >. In contrast, the higher 
velocity dispersion halos (which do not appear to undergo 
the ROI) show only mild initial increases in 2T r /T c f )1 due 
to their larger tangential velocities. For every halo, these 
increases correspond directly to the time when each halo 
deviates noticeably from its original spherical shape. 
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Fig. 6. — The ratio of the square of the velocity dispersion to the square of the circular velocity with radius, for each simulation at various 
times. Simulations represented are no dispersion (black solid line), la dispersion (red dotted line), 2a dispersion (green long— dashed line), 
and the 3a dispersion (purple dashed line). The timestep is labelled in each panel. Instabilities grow when the velocity dispersion is low, 
and stabilize when the velocity dispersion becomes comparable to the circular velocity. 



Following the peak in 2T r /T^, the onset of the ROI in 
the low-er halos produces a bar which torques the halo 
particles. These torques produce larger tangential veloc- 
ity dispersion, increasing < a^ > and thus decreasing 
2T r /T ( j ) to its final value of ~ 1.2. Note that all halos 
converge to this final value, regardless of the presence of 
the ROI. We plot the time to 5 Gyr to better portray 
the evolution of the halos; after this time, the values of 
2T r /T ( j > are fairly constant throughout the remainder of 
the simulation. 

Given the absence of the ROI triaxiality in the 2a and 
3cr halos, we can use the range between the peak values 
of the la and 2a runs to estimate the global anisotropy 
threshold for the ROI in our simulations. This estimate 
is a fairly loose one, since the la halo is not exactly 
spherical, nor in equilibrium. Nonetheless, the resulting 
value is 2.0 < 2T r /T^, < 2.97 for our halos, the lower 
range of which is consistent with previous estimates for 
spherical equilibrium models. 



The changes in anisotropy described above are not di- 
rectly associated with virialization. The lower panel of 
Figure [7] shows the quantity — 2T/W vs. time, where T 
and W are the sums of the kinetic and potential energies 
of each particle in a halo at each timestep. This quantity 
tracks the virialization process of the halo, which hap- 
pens on a different timescale than the anisotropy evolu- 
tion. Initially, each halo is a diffuse sphere expanding 
with the Hubble flow, so there is significant kinetic en- 
ergy and modest potential energy. As the simulation 
proceeds, the particles reach turnaround and begin to 
collapse, at which point the total kinetic energy reaches 
a minimum along with the quantity \2T/W\. The ROI 
then proceeds directly after the primary collapse of the 
halo, as can be seen by the coeval increase in 2T r jT^ 
(isotropization) with the rise of \2T/W\ (virialization) 
for each halo. However, isotropization due to the ROI 
continues after virialization, as can be seen in the broad 
peak in 2T r /T,p, which continues to evolve after \2T/W\ 
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Fig. 7. — (Upper Panel) Global anisotropy of each halo versus 
time, showing no dispersion (black solid line), la dispersion (red 
dotted line), 2a dispersion (green long-dashed line), 3<r dispersion 
(purple dashed line), and the 3a radial only dispersion (blue dot- 
dashed line). The ROI does not occur in systems with isotropic a 
of 2 or higher. (Lower Panel) Ratio of twice the kinetic energy of 
the halo to the potential energy versus time. The ROI occurs on a 
different timescale than the virialization process. 

has reached its final value. These different timescales 
emphasize the distinctness of these two processes. 

3.4. Anisotropy Profile 

Closely related to the global anisotropy is the 
anisotropy parameter, (3 (Eqn 1). The parameter /3(r) 
is also defined for a spherical system, but we use it here 
as an approximate indication of the amount of radial mo- 
tion in the orbits of our halos. When (3 = 0, the velocity 
dispersions <r r and are comparable and indicate an 
isotropic system. When /3 = 1, o~ r is much greater than 
cr^; the orbits in such regions are dominated by infall, 
"shell-crossing" , and in general highly radial orbits. We 
examine f3 calculated in radial shells to characterize the 
typical orbits of particles in a given region. The resulting 
anisotropy profile (3(r) can be combined with the shape 
information to isolate the dynamical effects taking place. 
We calculate (3{r) for each timestep using the values of 
ay and described in i )3.2l 

Our halos all settle into a common final anisotropy pro- 
file, characterized by a roughly isotropic core surrounded 
by an anisotropic, radially-dominated outer region, as 
shown in Figure [8j and as not ed in several other papers 
including Bar nes et alJ (|2005f ). However, we find that 
this profile occurs regardless of initial velocity dispersion, 
and whether or not the ROI has taken place. We can 
explain this behavior as follows: the anisotropy of the 
outer regions is expected, since material at large radii 
predominately consists of inward and outward moving 
particles that began with a high gravitational potential, 
hence having large characteristic radial velocities. In the 
center, velocities are isotropic for one of two reasons: ei- 
ther the ROI has occurred, transferring radial motions to 
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Fig. 8. — Anisotropy profiles for halos with no velocity dispersion 
(black solid line), la dispersion (red dotted line), 2a dispersion 
(green long-dashed line), 3<r dispersion (purple dashed line), and 
the 3a radial only dispersion (blue dot-dashed line). The horizontal 
black dotted line marks /3 = 1 for reference, and the vertical black 
dashed line marks the softening length. All halos develop isotropic 
cores and radially— dominated outer regions, regardless of initial 
velocity dispersions. 



tangential, or the motions are dominated by the initial 
isotropic velocity dispersion imposed by the initial con- 
ditions. In the latter case, any increase in radial velocity 
dispersion due to infalling shells is counteracted by the 
increase in tangential velocities as particles attempt to 
conserve angular momentum in response to being pulled 
inward by the infalling material. Therefore, both the 
large and small radius behavior suggest that the behav- 
ior of /3{r) increasing from to 1 should be quite typical, 
as we indeed find. 

Although the overall shape of (3(r) = — » 1 is similar 
among all halos, we do see a slight signature of the ROI 
in the steepness with which (3{r) increases. The higher- 
a halos have noticeably steeper transitions from f3(r) = 
— > 1, which could possibly be a result of the initial 
velocity dispersion dominating the core isotropy rather 
than the ROI. 

It is important to point out that a triaxial halo core 
cannot be purely isotropic, though it appears so in Fig- 
ure [5J This is due to our use of spherical radial bins 
to analyze a non-spherical system, and the fact that the 
standard-resolution halos do not have enough particles in 
the central region to properly resolve the orbital motions. 
An examination the Cartesian velocity dispersion tensor 
of the central regions confirms that the core is nearly but 
not exactly isotropic. Our use of (3{r) as a measure of 
anisotropy is still very useful, however, as it effectively 
tracks large-scale changes in anisotropy with radius, as 
can be seen from the figures. 

We track the time evolution of the anisotropy profile 
in Figure [5] for the high-resolution Oct halo. The char- 
acteristic shape of the profile develops very early on as 
particles at small radii collapse and undergo mixing. At 
the first timestep shown, the ROI has just begun, causing 
the innermost particles to enter box orbits, while the rest 
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Fig. 9. — Evolution of the anisotropy profile with time for the 
high-resolution, Oct halo. Red dotted line is t = 0.605 Gyr, green 
dashed line is t = 1.56 Gyr, blue dot-dashed line is t = 3.21 Gyr, 
black solid line is t = 12.8 Gyr. The black dotted line marks j3 = 1 
for reference. Each profile is normalized to the virial radius at 
that timestep. The characteristic shape of f3(r) forms early and is 
maintained throughout the halo's evolution. 



of the halo is still dominated by radial infall (red dotted 
line) . The near- isotropic core develops quickly and grows 
as more particles reach the center; by 1.6 Gyr, the ROI 
has operated on much of the center of the halo (green 
dashed line). The outermost regions consist of radially 
infalling material throughout the evolution of the halo. 
Note that the high-resolution halo has enough particles 
to resolve the orbital motions at the core. The /3(r) pro- 
file tends toward zero at small radii but an anisotropic 
component remains, due to the triaxiality of the halo. 

The time evolution of the anisotropy profiles of the 
higher-cr halo is comparable to those of the "cold" halo, 
in that the cores of the halos become isotropic at early 
times. However, while these halos appear to develop sim- 
ilarly, the core isotropy comes about differently for each 
halo; the higher-cr halos have pre-existing isotropic or- 
bits in the central region, preventing instabilities during 
collapse and maintaining a central isotropic core. In con- 
trast, the Otr halo develops its isotropic core via the ROI. 
Therefore, the characteristic shape of the anisotropy pro- 
file is independent of the shape of the halo, and can be 
formed either when there is an existing isotropic velocity 
dispersion, or via such a mechanism as the ROI if the 
halo is dynamically cold. 

This finding explains the results of iLu et al.l ((2006), 
who find that the universal 7^—1 inner halo density 
slope is formed for a variety of initial velocity dispersion 
distributions (i.e. radial, isotropic), particularly during 
an early "fast accretion" phase. Our Gaussian initial 
conditions are such that ~ 90% of the halo mass has a 
freefall time within 50% of the median. Thus our halos 
are lik ely to fall in the "fast accreti on" regime of | Lu et al.l 
(2006). Note, however, that unlike lLu etai] (|2006h we do 
not need to posit an additional source of isotropization 
for a core to be produced. The necessary isotropization 
occurs naturally either through the ROI or the conserva- 
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Fig. 10. — Fully evolved density and phase— space density profiles 
for halos with no velocity dispersion (black solid line), Ict disper- 
sion (red dotted line), 2ct dispersion (green long-dashed line), 3ct 
dispersion (purple dashed line), and the 3ct radial only dispersion 
(blue dot-dashed line). The vertical dashed line marks the radius 
below which softening and resolution begin to affect the profiles. 
All of the profiles are very similar, except that halos with higher 
velocity dispersion have shallower cores. 

tion of angular momentum. 

3.5. Density and Phase-Space Density Profiles 

Halos with various velocity dispersions show similar 
structural properties, whether or not they undergo the 
ROI. FigureQIIlplots the fully evolved density and phase- 
space density profiles of each standard-resolution halo. 
There is a noticeably shallower core in both the density 
and phase-space density profiles with higher a. Th is ef- 
fect was first reported bv lMerritt fc Aguilari (|1985[ ) and 
is not surprising, considering the higher initial angular 
momenta characteristic of systems with large o^. It is 
also important to recall from our tests in §2 that these 
medium-resolution simulations may not completely cap- 
ture the properties of the halo core. So, while these re- 
sults certainly make intuitive sense, high-resolution sim- 
ulations are needed for confirmation within 0.05r2oo- 

The evolution of the density and phase-space density 
profiles with time are shown in Figure [IT] for the high- 
resolution, zero dispersion halo. Each profile is normal- 
ized to the virial radius at each corresponding timestep. 
The basic shape of each profile is formed after 1.5 Gyr 
of evolution and grows outward with radius as the halo 
collapses. The outermost radii have yet to collapse at 
the present time. The final density profile (top panel) is 
fit well out to 10r2oo by an Einasto profile with a = 0.2. 
The NFW profile with c = 8.0 gives a poorer fit at the 
outer edges (see Fig. 1, top panels). 

Our halos exhibit a power-law in the phase-space den- 
sity (represented by p/c 3 ) with a slope of -1.7 that is 
relatively constant with time. The final halo power-law 
persists over ~ 2.5 decades in radius. The core phase- 
space density begins to develop its power-law structure 
at early times as well (Figure [TTJ bottom panel) . These 
results are consistent with other simulations and do not 
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Fig. 11. — Evolution of the density and phase— space density 
profiles with time for the high-resolution, Oct halo. Red dotted line 
is t = 0.605 Gyr, green dashed line is t = 1.56 Gyr, blue dot-dashed 
line is f = 3.21 Gyr, black solid line is t = 12.8 Gyr. All density 
curves are normalized to P200 at t = 12.8 Gyr to better show the 
density profile evolution. The profile shapes are set early on and 
evolve radially outwards with time. 

seem to be affected by the ROI. 

3.6. Scale Lengths 

iBarnes et al.l (|2005|) developed a semi-analytic model 
of collapsing halos which includes a physical represen- 
tation of the radial orbit instability. This model varies 
the initial perturbation velocities of the collapsing halo 
in such a way that an isotropic core is produced, with 
orbits becoming more anisotropic with increasing radius. 
This model produces realistic NFW density profiles and 
power-law phase-space density profiles for a series of ha- 
los. In addition, the model shows a direct correlation 
between the anisotropy radius (r a ), defined as the radius 
where (3 = 0.5, and the density scale length (r s , here 
defined as the radius where 7 = -2). This correlation 
is assumed to be a result of the ROI acting to flatten 
the density cusp by torquing the orbits of the inner halo 
particles and forming an isotropic core. 

While we showed in §3.3 that the characteristic 
anisotropy profile is not always due to the ROI, we still 
expect a correlation between r a and r s . Within r a , tan- 
gential velocities are large, and the resulting high angu- 
lar momenta should keep particles from penetrating the 
very center of the halo. This angular momentum barrier 
should lead to a flatter core whose size r s is physically 
related to the size of the isotropic region. 

To test the conclusion presented in IBarnes et al.1 
(|2005f ). in Figure [T^] we plot the relation between r a 
and r s as a function of time for the standard-resolution, 
zero-dispersion halo. The numbers on the figure indi- 
cate the time in Gyr at that particular point. While 
our N-body simulations do not confirm that the radii 
are equal, the two scale lengths are clearly proportional 
to each other. While there is considerable scatter, the 
points gravitate toward a line with slope ~ 2.5, suggest- 
ing that the anisotropy radius and the scale radius have 
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Fig. 12. — Anisotropy radius (r a /r20o) versus density scale radius 
{T s /r20o) for the standard resolution, Oct halo. These radii appear 
to be cor related with a slope of 2.5 (red line). The semi-analytic 
model of Barnes et al. (2005) gives a slope of ~ 1 (dotted line). 
Numbers on the plot indicate the time in Gyr at that particular 
point. 



some connection in their development. One reason for 
the difference in th e slope of our relation versus that of 
IBarnes et al.l (|2005f ) could be their assumption that the 
ROI does not change the energies of the particles, only 
the angular momenta. Our N-Body simulations allow for 
both of these quantities to change, allowing for a more 
accurate portrayal of the ROI. Thus the possibility re- 
mains that the ROI is the mechanism that produces the 
characteristic anisotropy profile for a "cold" collapsing 
system. 

4. CONCLUSIONS 

We have conducted N-body simulations of a variety of 
isolated, collapsing dark matter halos to investigate the 
role of the radial orbit instability on the final halo struc- 
ture. The ROI plays a role in determining the eventual 
shape, density profile, and phase-space density profile of 
a halo. The ROI is suppressed in halos with isotropic 
initial velocity dispersions, but halos with purely radial 
velocity dispersions appear to undergo the instability. 

Many numerical simulations have produced triaxial 
dark matter halos dFrenk et al.llT 988; Wa rren et al.lll992l 
iThomas et "ail 119981: Uing fc Sutol I2002D . Observational 
evidence for triaxial halos is also app earing, from low sur- 
face brightness and dwarf galaxies dBureau et al.l 119991 : 
ISimon et al.l 120051: iHavashi et al.ll2006l) and f rom galaxy 
clusters (|Oguri et all 120031 : iLee fc Sutdl2004r ). Such evi- 
dence for the existence of triaxial halos strengthens our 
hypothesis that dynamically cold halos undergo the ROI 
during collapse and therefore obtain a triaxial shape. A 
warmer system would not become unstable and would 
have a final shape that is too spherical to explain the 
observed structures. While the linear overdensities in 
the early universe ar e themselves triaxial, not spherical 
(|Bardeen et al.|[l986f ). and hence also contribute to the 
triaxiality of the resulting halos, the ROI significantly 
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affects the inner, observable parts of the halo. 

The characteristic anisotropy profile seen in all of our 
halos, and in others in the literature, is independent of 
the ROI, and appears to be a result of the general process 
of halo collapse. However, in cases where the ROI does 
occur, it is likely to be the cause of such profiles. In these 
cases, there is a link between the density scale radius 
and the anisotropy radius, indicating a causal relation 
between the two. Further semi-analytic and N-Body 
simulations may help to clarify this potential connection. 

While simulations of isolated collapsing halos are use- 
ful for trying to understand the effects of the ROI, it is 
well known that in realistic cosmological settings, halos 
form hierarchically. Their evolution is marked by peri- 
ods of gentle accretion, typically in the form of minor 
mergers, occasionally punctuated by major mergers. To 
understand the role and relevance of the ROI in halos 
evolving in this way, we are presently in the process of 
carrying out and analyzing simulations of h alos subject 
to co ntrolled major and minor mergers (cf. (|Poole et al.l 
2007)). At this early stage, we are finding that the com- 
plex dynamics at play during major mergers makes it 
difficult to ascertain straightforwardly whether the ROI 



plays a significant role during the subsequent relaxation 
process; a priori, we would speculate that it does not. 
However, there are indications that the ROI is opera- 
tional during periods of quiescent accretion, by acting 
upon nearly radial tidal streams associated with disrupt- 
ing subhalos, which become isotropic. We speculate that 
the weak tides induced by these weak structures in fact 
both seeds and reinforces the instability. 
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